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ABSTRACT 



Observational data and theoretical calculations show that significant small-scale substructures are present in dark molecular clouds. 
These inhomogeneities can provide precious hints on the physical conditions inside the clouds, but can also severely bias extinction 
measurements. We present Nicest, a novel method to account and correct for inhomogeneities in molecular cloud extinction studies. 
The method, tested against numerical simulations, removes almost completely the biases introduced by sub-pixel structures and by the 
contamination of foreground stars. We applied Nicest to 2MASS data of the Pipe molecular complex. The map thereby obtained shows 
significantly higher (up to 0.41 mag in Ak) extinction peaks than the standard Nicer (Lombardi et al. 2001) map. This first application 
confirms that substructures in nearby molecular clouds, if not accounted for, can significantly bias extinction measurements in regions 
with Ak > 1 mag; the effect, moreover, is expected to increase in more distant molecular cloud, because of the poorer physical 
resolution achievable. 
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1. Introduction 

Although the details of star and planet formation are very lit- 
tle known, there is a large consensus that these objects are 
created inside dark molecular clouds from the contraction and 
fragmentation of dense, cold cores. As a result, in the last 
decades molecular clouds have been studied in detail using many 
different te chniques, from optical number counts (Woll|[ 1923t 
iBokl Il937h , to radio observations of carbon monoxide (CO) 
molecules dWilson et alj 119701: iFrerking et all Il982h. to near- 
infrared (NIR) extinction measurements (^ Lada et al Il994 . 

Molecular hydrogen and helium represent approximately 
99% of the mass of a cloud, but the absence of a dipole moment 
in these molecules makes them virtually undetectable at the low 
temperatures (T ~ 10 K) present in these objects. Hence, the 
(projected) mass distribution of dark clouds is usually inferred 
from the distribution of relatively rare tracers (such as CO or 
dust grains) under the assumption that these are uniformly dis- 
tributed in the cloud. 

As pointed out bv lLada et alj (1 19941) . the reddening of back- 
ground stars in NIR bands is a simple and reliable method to 
study the dust distribution and thus the hydrogen column den- 
sity. Compared to optical bands, NIR bands are less affected by 
extinction and are less sensitive to the physical properties of the 
dust grains (Mathis 1990), and thus their color excess can be 
directly translated into a hydrogen column density. In a series 
of papers we reconsidered and improved the NIR color excess 
(Nice) technique, by ge neralizin g it to use three or more bands 
(Nicer, see Lombardi & Alves 2001 ) and by taking a maximum- 
likelihood approach (,Lombardi.2005i) . 
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Although the Nice and Nicer techniques hav e been very 
successful in studying molecular clouds (see, e.g., ' Alves et al] 
2001), they are plagued by two complications: small-scale inho- 
mogeneities and contamination by foreground stars. NIR stud- 
ies typically assume a uniform extinction over (small) patches 
of the sky: however, in presence of substructures such as steep 
gradients, unresolved filaments, or turbulen ce-induced inhomo - 
geneities (iLada etaI.lll994tlLarso n 1981: Hever & Bruntll2004 . 
the background stars used to estimate the cloud extinction would 
not be uniformly distributed in the patch, but would be preferen- 
tially observed in low density regions. Similarly, in presence of 
foreground stars, their relative fraction increases with the dust 
column density. Unfortunately, both effects will bias all color 
excess estimators toward low column densities: moreover, the 
bias will be especially important in the dense regions of molec- 
ular clouds, which is where the star formation takes place. 

In this paper we describe in details the repercussions of small 
scale inhomogeneities and foreground stars in extinction mea- 
surements, and propose a simple method to obtain estimates of 
the column density in presence of both effects that, under sim- 
ple working hypotheses, is asymptotically unbiased. Our method 
does not rely on any (usually poorly verified) assumption regard- 
ing the small scale structure of the cloud: moreover, it can be ap- 
plied to a general class of NIR extinction measurements (marked 
spatial point processes). 

The paper is organized as follows. In Sect. |2] we introduce 
a general smoothing technique virtually used in all cases to cre- 
ate smooth, continuous maps from the discrete pencil-beam ex- 
tinction measurements caiTied out for each background star: we 
also show that in two typical applications (moving average and 
nearest neighbours smoothing) a bias is expected. In Sect. |3] we 
propose a new method that can correct this bias without making 
any assumption on the underlying form of the cloud substruc- 
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ture. This new method has been tested with numerical simula- 
tions, as described in Sect. H) and with a real-case application, 
presented in Sect. |5] We discuss and comment the results ob- 
tained in Sect. |6l and finally we briefly present our conclusions 
in Sect, in 



2. Cloud substructure 

Measurements of the reddening of stars observed through a 
molecular cloud provide estimates of the cloud column densi- 
ties on pencil beams characterized by a diameter of the order of 
a fraction of milliarcsec. However, these high resolution mea- 
surements are strongly undersampled and are affected by large 
uncertainties due to the photometric errors and to the intrinsic 
scatter of star colors in the NIR. Smooth and continuous extinc- 
tion maps are normally obtained by interpolating and binning 
the various pencil beams. Di fferent authors u se dififerent recipes 
to smooth the data (see, e.g.. lLombardil20d 2^ for a description of 
several interpolation techniques), but in most cases the interpola- 
tion follows a standard scheme. Consider /T-band extinctiorQ 
measurements {A„] obt ained, for example , fr om the color excess 
of background (see iLada et all 1 19941 or iLombardi & AlvesI 
12001 ). The extinction A at any location in the sky is evaluated 
using a weighted average 



(1) 



The weights {w„] are usually chosen to be significantly differ- 
ent from zero only for stars angularly close to the given location 
of the map (for example, Cambresy et al. ,2002 assign a unity 
weight to the nearest neighbors). Clearly, Eq. ([T) does not 
include slightly more complex situations where, for example, 
one takes the median of th e extinctions nieasur ed from objects 
nearby a given position (e.g. lDobashi et al.l2008) : however, most 
of the discussion carried out for the weighted average actually 
applies also to median or related estimators. 

The binning in Eq. ([T) washes out the cloud substructure on 
scales smaller than the typical size where the {w,,] are signifi- 
cantly different from zero, but this is needed in order to have 
smooth maps and to increase the signal-to-noise ratio. However, 
Eq. ^ also introduces a significant bias on the estimated col- 
umn density. Suppose that, in the region of the sky that we are 
investigating (i.e., in the area covered by the stars used to es- 
timate A) the column density has significant variations. Because 
of this, the local density p{x) of stars is not homogeneous though 
the cloud, but rather follows the scheme 



p{x) =polO 



-ak\A(x) 



(2) 



where po is the density of stars where no extinction is present, a 
is the slope of the number counts, and k,i = A ^/A is the extinction 
law in the band A considered. This effect, which is at the origin of 
the historical number count method to measure column densities 
in molecular clouds (Wolf 1923; Bok 1937), also induces a bias 
in Eq. ([TJ, since regions with smaller extinctions and thus higher 
density of background stars will, on average, contribute more to 
the sum in Eq. ([T]i than regions with large extinctions. Note that 
since the color excess method requires measurements in at least 
two different bands, one should replace Eq. ^ with a version 



' For simplicity, in this paper we drop the subscript K normally used 
in literature to denote the A'-band extinction Ax', it is also obvious that 
the same method applies to extinction measurements referred to any 
band (e.g., the visual extinction Av). 



that provides the expected density of stars observed in all bands 
required for the application of the method (e.g., H and K). For 
equally deep observations on all bands, the result has the same 
form of Eq. (|2|i, where k,\ is the extinction law of the bluer band 
considered (e.g., H). 

The main issue considered in this paper is the bias intro- 
duced by unresolved substructures in Eq. ([U. In general, the bias 
should be defined here as the difference between the expected, 
mean value of (A) and the true column density A at the same po- 
sition. However, clearly every smoothing technique introduces a 
bias because of the smoothing itself, and this bias is normally ac- 
ceptable if it does not introduces a systematic effect on the total 
column density. In other words, a required property is 



J A(x)dx = I^J A(x)dx^ = J {A(x))dx 



(3) 



Clearly, in order for Eq. Q to hold, the difference (A) - A must 
be of alternating sign. For our purposes, thus, it is more useful 
to define the bias as the difference between the expected mean 
value (A) and the same quantity that we would theoretically ex- 
pect in presence of a uniform density of stars, i.e. by ignoring 
the dependence of the density of stars on the local extinction de- 
scribed by Eq. (|2]l [see also below Eq. O]. This definition is 
useful, since it isolates the standard effects of a smoothing from 
the effects introduced by the non-uniform sampling of the cloud 
structure. In addition, a column density estimator A(jc) that is un- 
biased according to this definition will also be unbiased accord- 
ing to Eq. (O (provided the weights w„ are spatially invariant). 
For, when the density of background stars is uniform, {A{x)) can 
always be written as a convolution of the true column density 
map A (x) with some kernel K normalized to unity (iLombardil 
I2002h . 

In the foUowing, we will calculate explicitly this bias in two 
common smoothing schemes. 

2.1. Moving average smoothing 

In this section we will consider a simple weighting scheme in 
Eq. ([T) where each weight w„ = w(x„; x) is a simple function of 
the location x„ of the «-th star (typically, w{x„;x„) will depend 
only on the modulus \x„ - x\). 

As discussed above, our aim is to compare the expected col- 
umn density estimate with the one that we would obtain in ab- 
sence of any selection effect in the number of background stars. 
We need thus to evaluate two averages: (i) {A{x)), where A is cal- 
culated according to Eq. ([T]i and the ensemble average is taken 
over all positions of stars with the density given by Eq. ^ and 
over all individual column density measurements A„; and (ii) 
the average A(x), which is basically the same quantity evaluated 
with a constant density of background stars po. 

In principle, both average s can be calculated analytically us- 
ing the method described by iLombardi & Schneider! (1200 11: see 
below Sect. [3]i; in practice, for the purposes of this section a 
simple approximation can be used provided that in the weighted 
average of Eq. ^ a relatively large number of column densities 
are used with weights significantly different from 0. In this case 
we find for the first average 



<A(^)> 



j A{x')w{x';x)p(x') dx' 
j w{x';x)p(x')dx' 



(4) 



As shown by Eq. p{x') is a simple function of A{x'). This 
suggests that the equation can be recast in a simpler form by 
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defining pa(A;x), the (weighted) probability that a point with 
extinction A is used to evaluate A{x), the column density at x: 



PAiA;x) = 



/ w(x';x)6{A-A(x'))dx' 
J w(x';x) dx' 

By using Eq. ^ in Eq. (HJl we find then 



(Mx)) = 



/ Apa(A; x)p(A) dA J Apa{A; a:)10-«*^"^ dA 
J PA{A-x)p{A)dA ~ //?A(A;jc)10-«*i^dA 



(5) 



(6) 



In presence of a uniform distribution of stars, p(x) - po, we 
would instead obtain simply 



Aix)^ 



J A(x')w(x' ; x) dx' 
J w(x';x) dx' 



dx' r 

— = J ApA(A-x)dA. 



(7) 



The two values presented in Eqs. (|6]l and O do not differ signifi- 
cantly provided the scatter of A in the patch considered around x 
is much smaller than (akA In 10)"'. For example, if we consider 
the H band of a typical region close to the Galactic plane with 
a standard reddening law (Rieke & Lebofskv 1985), we have 
a ^ 0.34 mag ' and kn - 1.55 dlndebetouw et al.ll2005b . so that 
the maximum scatter allowed in Ak to have a negligible bias 
is <s 0.82 mag. Hence, clearly we need to be concerned by this 
bias in regions from moderate to large extinction. 

It is also interesting to evaluate the bias {A{x)) -A(x) in pres- 
ence of small variations of A within the region considered. In this 
case, the probability distribution pa{A;x) is peaked around the 
mean value of Eq. (|7]), and we can expand to second order the 
two exponential functions present in the numerator and denomi- 
nator of Eq. (|6]l. After a few manipulations we obtain then 

{A(x)) - A(x) ^ - aki In 10 J pa{A; x)[A - A{x)f dA 
^ w{x';x)l:i^(x';x) dx' 



J wix';x) dx' 



(8) 



where /3 = ak^ln 10 and A{x';x) = A(x') - A(x). Hence, the 
difference between the two values is proportional to a weighted 
average of the scatter of A around its mean value A(x) at x, a 
quantity that, as we will see below, can be directly estimated 
from the data. Note that the bias is, to first order, quadratic on 
A(x'; x) and thus will be particularly severe when steep gradients 
are present in the underlying column density A{x). 

2.2. Nearest neighbor(s) 



ICambresv et al.l (|2002|) suggest to use a different prescription to 
make smooth extinction map from the individual column density 
measurements. For each point in the map, they take the average 
extinction of the angularl y closest stars (nearest neig hbours 
interpolation). As argued bv 'Cambresv et al] (I2005L 1200 6). this 
method can potentially alleviate the bias introduced by the vary- 
ing background density of stars described by Eq. (|2]l, because 
measurements from low density regions will mostly appear iso- 
lated and will thus be used for relatively large patches of the sky. 

Interestingly, using the technique described by iLombardl 

( 12002 ) it is possible to obtain an analytical estimate for the aver- 
age of A in the smoothing scheme proposed: 



{A{x)) = A(x')K(x'; x)p{x') dx' 



(9) 




Fig. 1. The bias of the nearest neighbors interpolation on a toy 
model where the true extinction is A(x) -O.l mag for jci < and 
0.9 mag for xi > 0. The solid lines show the average measured 
extinction of the nearest neighbours method, calculated using 
Eq. (|9]l for various number of neighbours A^. The dashed lines 
show the extinction that one would measure if the density of 
stars were uniform on the whole field. Note that the solid lines 
are always below the dashed lines except at the extremes of this 
plot (where all estimators agree with the true value of A there). 
For this plot we set a - 0.34, k = 1 .55, and po - I. 



where the linear kernel K(x';x) is given by 



K(x';x) = y Itl^Ll^ 



N 



(10) 



In this equation, is the number of stars used in the nearest 
neighbors interpolation and fi{x';x) is the average number of 
stars observed in B{x'; x), the disk centered in x and of angular 
radius |jc - jc'l: 



Hix'-x)= f pix")dx" 

JB{x';x) 



(11) 



As shown by iLombardil (l2002h . the kernel K(x',x) is normal- 
ized, i.e. 



K(x';x)p(x')dx' = 1 , 



(12) 



an obvious result if one consider the measurement of the extinc- 
tion on a uniform cloud where A(x) - const. In our case, this 
property can also be proved directly from Eqs. (fTOl i and ( fTTT i. 
First, note that fi(x';x) - p.{r\x) depends only on x and on 
r - \x - x'\, and thus the same applies to K{x';x) - K{r,x). 
This suggests that we can recast the integral of Eq. ( fT2l l as 



f 

Jo 



dfiinx) 

K(r, x) — dr 

or 

N-l 



-y- r 



-'^('-'[M'-;jc)r^4^dr=l. (13) 

or 



The fact that K{x';x) only depends on x and r - \x-x'\ (and 
not on the direction of the vector x - x') also implies that the 
nearest neighbors interpolation suffers from the bias discussed 
in this paper. Indeed, in the integral of Eq. (|9]) the kernel K gives 
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Fig. 2. Similarly to Fig.[T] but for the median-nearest neighbors 
interpolation. Again, the solid lines, representing the expected 
mean measurements, are always below the dashed lines, repre- 
senting the same measurements in the ideal case where the den- 
sity of background stars is uniform within the whole field. Note 
that because of the properties of the median, for large values of 
the "transition" from A = 0. 1 mag to A = 0.9 mag is faster 
than in Fig.[T] 



the same "weight" to all points on circles centered in x, irrespec- 
tive of the specific value of p on the various points of the circles 
(what matters here is only the average value of p on a circle, 
and not the specific value on the various points). Hence, we do 
not expect that this technique can solve the bias introduced by 
the correlation between the extinction A and the density of back- 
ground stars p expressed in Eq. (|2]i, especially if steep gradients 
are present in the intrinsic cloud column density. We also ex- 
pect that the bias present in the nearest neighbours interpolation 
increases with the number of neighbours A^. Indeed, when in- 
creases, the "size" of the kernel K(r, x), i.e. the range in r where 
the kernel is significantly different from 0, also increases because 
of the effects of the polynomial terms in Eq. ( fTOb . This, as shown 
by Eq. (|9]), implies that the kernel averages out more distant re- 
gions in the sky, where differences among the local densities can 
be significant. 

In order to better understand the points mentioned in the pre- 
vious paragraph, it is useful to look at simple example. Consider 
a cloud characterized by a Heaviside column density: 



A{x) 



[Ai if jci < , 
1 A2 if xi > . 



(14) 



In this case we can basically carry out all the calculations analyt- 
ically, and obtain for each value of xi the expected average mea- 
sured extinction (A). We do not report here the relevant equa- 
tions, but rather show the results obtained in a typical case in 
Fig-El As shown by this plot, the nearest neighbours are clearly 
biased: for example, the various curves are not symmetric around 
xi - 0, and in addition 

AilO''*^^' +A2lO"'*'^^ A1+A2 
(A) = r-, r-. < ■ (15) 

In order to better understand the bias, we also plot in Fig. [T] the 
expected results in case of a uniform density, obtained assuming 
a = in Eq. (|2]) (dashed lines): in this case the curves are per- 
fectly symmetric around xi = and (A) = (Ai + A2)/2. A com- 
parison of the solid and dashed lines also confirms that curves 



Fig. 3. Similarly to Fig. [T] but for the moving average smooth- 
ing with a Gaussian window function with different dispersion 
parameters cr. As for Figs. [T] and |2l the solid lines, represent- 
ing the expected mean measurements, are always below the 
dashed lines, representing the same measurements in the ideal 
case where the density of background stars is uniform within the 
whole field. 



with large A^, being less steep, are biased on a larger interval 
around xq = 0. This is expected, since as shown by Eq. ( fTOt the 
intrinsic size of the smoothing kernel K increases with A^. 

One might wanders if the bias decreases by using a median 
estimator ins tead of a simple mean, as done by ICambresv et al.l 
(|2005L|2006|) . These authors still identify, for each sky position 
X, the A^ nearest neighbour stars, and then evaluate a simple me- 
dian of the relative extinction measurements. Luckily, it is possi- 
ble to eval uate the exa ct statistical properties of the median esti- 
mator (see lLombardil|20 02 , Appendix A). For this purpose, it is 
useful to evaluate the probability distribution pn(A;x) that one 
of the A^ neighbours around the position x has a column density 
A: 

Pm(A;x)^ fdije ''^ [ dx' p(x')6{A - A(x)) , (16) 

where B(fi; x) is the disk centered at x and with radius r defined 
by yU = yu(r; x) (note that since //(r; x), for x fixed, is a mono- 
tonic function of x' this definition is unique). The cumulative 
distribution associated with p^iA; x) is given by 



Pn{A-x)^ I dpt-^''^v{fi,A;x) 



(17) 



In this equation, v(yU, A; x) is the integral of p{x) carried out over 
all points of B(/i; x) where A{x') < A: 



v(jU,A; jt) 



p(jc')dx' 



(18) 



B(ii;x)r\{x'\A{x')<A] 



Note that v(p.,A\x) — > /i as A — > 00. Finally, the median of the 
N -2k-\ nearest neighbours is then provided by the expression 



Pm{A) 



kV^ \pi,{A-x)[Pn{A-x)T'[\-Pn{A-x)] 



k-\ 



(19) 



Although it is difficult to obtain a general expression for the 
bias introduced by the median-nearest neighbours estimator, it 
is clear that to first approximation the same arguments discussed 
above for the simple mean apply (note also that, trivially, when 
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N - I the two estimators are identical). However, for the simple 
example considered above, we can actually carry out the calcu- 
lations using Eqs. ( fT6HT9l ) and show that the bias is still signifi- 
cant. Figure |2] summarizes the results obtained in the toy model 
described by Eq. (fT4l i. and proves that no real improvement is 
obtained from the use of a median estimator (cf. Fig.[T]). In par- 
ticular, a severe bias is still present, and its amount increases 
with A^. For comparison, in Fig. [3] we also show the results ob- 
tained from the simple moving average smoothing discussed in 
Sect. 12. II interestingly, the plot is very similar to the one for a the 
simple average nearest neighbors interpolation (Fig.lTJ. This, in 
reality, is not too much surprising because the expected average 
values (A(x)) measured with these two techniques can be de- 
scribed in terms of a simple convolu tion with some a ppropriate 
kernels [cf. Eq. dUi and ([foli; see also Lombardill2002ll . 

In summary, a careful analysis of the nearest neighbors 
shows that this interpolation technique, for the specific case of 
extinction measurements, is still affected by a significant bias. 
At least for the case of a simple mean estimator, the origin of the 
problem can be traced to the different behaviour of the kernel 
K{x; x') with respect to the density p(x') in Eq. (|9]i: in particu- 
lar, the kernel K does not respond directly to variations of p, and 
instead depends only fx(x';x), i.e. on the total "mass" within a 
disk of radius r = (x - x'}. 



3. Nicest: over-weighting high column-density 
measurements 

As shown by Eq. ([8]l, the bias discussed here strongly depends 
on the small-scale structure of the molecular cloud through the 
probability distribution pAiA). Although on large scales many 
models predict a log-normal distribution for p^, clearly we 
should expect significant deviations from the theoretical ex- 
pectations on the small scales often investigated in NIR stud- 
ies. In addition, even at large scales the superposition of dif- 
ferent cloud complexes, or a significant "thickness" of a cloud 
along the line of sight, c an produce significant deviations fro m a 
log- normal distribution (|Vazquez-Semadeni & GarcialuOOlt see 
also lLombardi et al.l200 6lfor a case where the log-normal distri- 
bution is not a good approximation). Hence, we need to be able 
to correct for the substructure bias in a model independent way. 

3.1. Statistical properties of tiie moving average smootiiing 

In order to better understand, and then fix, the bias present in 
Eq. ([T]) in the case o f the moving average srnoothing, we use th e 
theory developed in lLombardi & Schneide3(l200llEooll2003h . 
The key equations to obtain the ensemble average of A{x) are 
summarized below: 



function does not effect Eq. ([T]i. For simplicity, in the fol- 
lowing we will assume, without loss of generality, that the 
weight function is normahzed according to 



w{x'\x)p{x')Ax' - 1 



(24) 



The quantity Q(s\x) can be shown to be related to the 
Laplace transform of p„{w;x), the probability distribution 
for the values of the weight function w(x'; x) with x fixed. 
P(x) is the probability that no single star is present within the 
support n^ix) of w{x'; x), defined as n„{x) = {x' \ w(x'; x) > 
0). Hence, P(x) can be simply evaluated from 



P(x) - exp 



p(x')dx' 



(25) 



Note that P(x) - for weight functions with infinite support 
(such as a Gaussian). 

C(w;x) is a simple Laplace transform of Y{s;x). This is the 
key quantity that enters the final result (l23T l together with 
the original smoothing function w(x'; x) and the star density 
pix'). 

A key parameter in the calculation of C{w), and thus on the 
final result {A{x)}, is the so-called weight number of objects 
Nix), defined as 



Nix) = 



[1 



Pix)] J [wix';x)fpix')dx' 



(26) 



where as usual wix';x) is taken to be normalized according 
to Eq. (l24l i. Informally, Nix) counts the number of stars that 
contribute (with a weight significantly different from zero) 
to the average Aix). 

In the limit Nix) » 1 it is possible to obtain a simple ex- 
pression for C(w; x), which for Pix) - takes the form 



Ciw;x) 



1 



[Nix)]- 



l + W (1 H- w)^ 



(27) 



This expression shows that to lowest order C - 1 - w + A/' ' , 
where both terms w and A/' ' introduce a relative correction 
of order AA"' in the final result {A{x)). 

The results summarized above rigorously prove the state- 
ments of Sect. 12.11 In particular, in the limit of a large weight 
number of objects Nix), we have C(w) ^ 1, and thus Eq. ( |23] l 
together with the normalization (l24l i reduces to Eq. (HJi. In prac- 
tice, numerical simulations show that Eq. ^ is already accurate 
for relati vely small weight numbers (o f the order of A/' ~ 10; see 
Fig. 3 of lLombardi & Schneidedl200lh . 



Qis;x)^ 




- l]pix')dx' 


(20) 


Yis;x) = 


exp[Qis; x)] 




(21) 


Ciw;x) = 


1 p 


e-'"y(s;jc)ds , 


(22) 


I -Pix) Jo 


<A(^)> = 


1 Aix')w(x' 


x)C{w{x';x),x)pix')dx' . 


(23) 



Let us briefly comment this set of equations: 

- The equations are invariant upon the transformation 
wix';x) i-> qwix';x), with q positive constant. This is ex- 
pected, since an overall multiplicative constant in the weight 



3.2. A fix for the inhomogeneities bias 

The bias of Eq. ([TJ is essentially due to a change in the density 
of extinction measurements as a function of A [Eq. ^]. In the 
framework of the moving average smoothing, it is thus reason- 
able to try to "compensate" this effect by including in the weights 
of the estimator ^ a factor 10"*"*: 



Wn - wiXn', X)10' 



ak lA.. 



(28) 



With this modification, we expect A to be unbiased provided that 
N is large (so that we can take C - 1) and that measurement er- 
rors on A„ can be neglected. Instead, in presence of significant 
errors on the individual column density estimations, we expected 
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Fig. 4. Similarly to Fig. [3] but for the Nicest moving average 
smoothing (see Sect. [3]). In this case a significant bias is observed 
only for cr = 0.5, which however correspond to a very small 
weight number N -2. 

A to be biased toward large extinctions. This happens because of 
the non-linearity of the introduced factor in A„: for example, A„ 
is symmetrically distributed around A„, 10"'^ ''*"A„ will be skewed 
against values larger than 10"*^ ''^"A,,. In summary, the modified 
weight of Eq. ( l28l l removes (to a large degree) the bias due to the 
cloud substructure, but introduces a new bias related to the scat- 
ter of the individual extinction measurements. While the former 
is basically unpredictable, the latter can be estimated accurately 
provided we know the expected errors on A„ . 

In order to better understand the properties of the estimator 
([TJ with the new weight, let us evaluate the average (A). For 
simplicity, initially we will ignore measurement errors and focus 
on the effects introduced by a non-uniform density. We first note 
that Eq. ( l28l l is equivalent to the use of a weight w„ — w'{Xn',x), 
where 



w'{x'\x) = w(x';x)/p(x') 



(29) 



We require the new weight w'{x'; x) to be normalized according 
to Eq. (|24] |. which implies 



w{x';x)dx' - 1 . 



(30) 



Interestingly, the normalization condition (l30l l for w does not de- 
pend on the density p anymore, and it is thus possible to choose 
for w a spatially invariant function w{x';x) - w(x' -x) (a typical 
choice would be a Gaussian, w oc e"'* ^-^l/^o""). We have then 

(Aix)) = J A(x')w'ix'- x)[l - w'(x'; x) + N'\x)]pix') dx' 
Aix')wix';x)[l-wix';x)/pix') + N-\x)]dx' . 



I' 



(31) 



After a few manipulations we can rewrite this expression in the 
form 



(Aix)) = A(x) - — r w\x'; jc)10°'*^'*<^''[A(jc') - A(x)] dx' 
Po J 

^A(jc) w^{x'\x)[A{x';x) 

Po J 



+ /3A^ix';x)]dx , (32) 



where in the second equality we have expanded the exponen- 
tial term to first order using the definitions following Eq. ([8]). In 
summary, the bias {A(jc)) - A{x) of the proposed estimator is 
composed of two terms: 

- A weighted average of A(jc'; jc). Since the weighted average 
uses w^{x';x) as weight, this term is not expected to vanish 
identically unless the weight is a top-hat function [cf. Eq. (|8]l, 
where instead the weighted average involves w(x';x) and 
therefore no linear contribution in A is present]. Still, we do 
not expect any systematic bias related to this term, because 
on average A(jc'; jc) by construction has alternating sign. 

- A weighted average of A^{x'; x). This term is manifestly neg- 
ative and therefore introduces a bias in the result. The bias is 
very similar to the original bias discussed in Eq. ([8]): it is 
proportional to /3 = ak,\lnlO, and depends on the scatter 
A{x';x) = A{x') - A{x) of the local column density with 
respect to its average value. However, there is a significant 
difference: the bias is inversely proportional to polO"'*''**'^*, 
i.e. to the local average density within the weight function; 
moreover, the average is taken over instead of w. These 
two differences, together, indicate that the bias of Eq. (|32] |. 
compared to the on of Eq. ([8]l, is smaller by a factor ~ N(x). 

In other words, the modification of the weight operated by 
Eq. ( |28] ) reduces the bias present in the original simple mov- 
ing weight average by a factor Nix), which in typical cases is 
significantly larger than unityQ On the other hand, the fact that 
the method is ineffective then N is low is evident from the def- 
inition of w„: in particular, as p — > we expect that only a 
single star contributes to the averages of Eq. ([T) with a weight 
w„ - w(jc„; x)10"*'^" significantly different from zero, but then 
trivially A = A„ and no correction is effectively performed^ 

In Fig.|4]we show the effects of the proposed smoothing tech- 
nique on the toy-model discussed in Sect. 12.21 A comparison 
with Fig. [3] shows that the Nicer interpolation is very eff'ective 
in reducing the bias of the simple moving average interpolation, 
especially for relatively large values of the smoothing factor cr 
(and thus of the weight number of objects AO- 

We now consider the effects of measurement errors in the 
modified estimator. For simplicity, we consider the limit of a 
large number of stars (N » 1), so that C(w) ^ 1 (in any case, 
as shown above, the technique proposed is inefficient when N is 
of order of unity). In this case, a simple calculation shows that 



Be 



(33) 



where {cr^} are the variances on {A„}. In other words, the method 
proposed in this paper introduces a small bias toward large ex- 
tinction values. In order to estimate the order of magnitude of 
Ben, we note that the median error in K band extinctions on 
Nicer column density estimates from the 2MASS catalog is 
typicall}|3 cr ^ 0.13 mag, and that less than 1% of stars have 
cr > 0.2 mag. Taking cr„ - 0.13 mag for all stars in Eq. d33T l, 



^ In reality, a detailed calculation shows that the key parameter here 
is Ntsix), which is defined an alogously t o N(x ) but with w replaced 
by Weff = wC(w). As shown in lLombardil ( l2002h . N^ff > 1, so the bias 
never increases. 

' In this example N while N^s — > 1 , which agrees with state- 
ment of the previous footnote. 

The expected error on A„ is calculated using the standard Nicer 
techniques, and thus includes both the typical 2MASS photometric un- 
certainties and the star color scatters as measured in control fields where 
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we find a bias (A) - A ^ 0.02 mag (always in the K band). 
Interestingly, although very small this bias, in contrast to the one 
described by Eq. ( |32] |. can be corrected to first order because its 
expression involves only known quantities. 

In summary, we propose to replace the estimator ([T) with 



(34) 



where w„ = w(jr'; jr)10'*'^". This new estimator, that we name 
Nicest, significantly reduces the bias due small-scale inhomo- 
geneities in the extinction map and is (to first order) unbiased 
with respect to uncertainties on A„. 

If one is ready to accept a small bias on A for extinction mea- 
surement uncertainties, it is also possible to use Eq. ( l34b without 
the last term. In this respect, we also note that the correcting 
factor present in Eq. ( l34b is independent of the local extinction, 
and only depends on the individual errors on the measurements 
(these typically are a combination of the photometric uncertain- 
ties and of the intrinsic scatter of colors of the stellar population 
considered). As a result, it is not unrealistic to expect that this 
factor is approximately constant within the field of observation, 
and that it is equally present in the control field used to fix the 
zero of the extinction. 

Finally, we note that in Nicest a central role is played by 
the slope of the number counts a, operationally defined through 
Eq. In reality, however, it clear that Eq. ^ cannot strictly 
hold, because there is an upper limit on the (un-reddened) mag- 
nitude of stars. Hence, depending on the limiting magnitude of 
the observations, there is an upper limit to the measurable extinc- 
tion of stars. As a result, in regions where the extinction is larger 
than this limit, the density of background stars vanishes, with the 
result that in these cases even Nicest cannot be unbiased. The ex- 
act value of the upper measurable extinction depends on many 
factors, including the waveband A and limiting magnitude m^ax 
of the observations, and the distance d of the molecular cloud. 
An approximate relation for the maximum extinction Amax is 



1 

kx 



-5 log 



10 



10 pc 



(35) 



where M„iax is the maximum absolute magnitude of stars in the 
band considered. The maximum extinction A^ax should be eval- 
uated in each specific case; for example, for t ypical 2MASS ob- 
serva tions, using m^ax - 15 mag, kn - 1.55 dlndebetouw et al.l 
120051), and M,„ax ^ -5 mag we obtain Aniax - 1 mag for a cloud 
located at 100 pc. 

3.3. Foreground stars 

Foreground stars do not contribute to the extinction signal, but 
do contribute to the noise of the estimators ([TJ and ( [34l i, and 
thus whenever possible they should be excluded from the analy- 
sis. Usually, foreground stars are easily identified in high {Ak > 
1 mag) column density regions, but are almost impossible to dis- 
tinguish in low and mid extinction regions. So far we assumed 
that all stars are background to the cloud, but clearly in real ob- 
servations we should expect a fraction / of foreground objects. 
In principle, if this fraction is known, we could correct the col- 
umn density estimate A into A/(l - /): in this way, we would 



the extinction is negligible. Hence, this error includes the effects of dif- 
ferent stellar populations on the scatter in the intrinsic star colors, as 
long as the various stellar populations are represented in the control 
field. 



obtain an unbiased estimator of the column density at the price 
of an increased noise (due to foreground objects). In practice, 
the fraction of foreground stars is itself an increasing function of 
the local column density, i.e. / - /(A), because of the decreased 
density of background stars in highly extincted regions. Hence, 
the simple scheme proposed above is not easily implemented. 

Interestingly, Nicest perfectly adapts to the presence of fore- 
ground stars. Suppose that in regions with negligible extinction 
a fraction /o = /(A = 0) < 1 of stars is foreground to the cloud: 
then, in our notation, we can rephrase this by assigning a non- 
vanishing probability to Pa{A - 0), i.e. by treating foreground 
stars as a special case of substructure (much like "holes" in the 
cloud). As we consider higher density regions, the fraction /(A) 
of foreground stars increases but, because of the correcting fac- 
tor in w„, we still expect to measure (1 - /o)A: in other words, 
the estimator A /(I - /o), with A given by Eq. ( [34] l is expected to 
be unbiased both for small-scale inhomogeneities and for fore- 
ground contamination. Note that coiTecting factor (1 - /o)"' is 
usually very close to 1 for nearby molecular clouds, and can be 
easily evaluated by comparing the density of foreground stars 
(easily measured in highly extincted regions) with the total den- 
sity of stars in regions free from extinction. 



4. Simulations 

In order to test the effects of substructures and the reliability of 
the method proposed here in a real case scenario we performed 
a series of simulatio ns. We started f rom a 2MASS/Nicer map 
of the Pipe nebula dLombardi et al.l 12006). a nearby complex 
of molecular clouds seen in projection to t he Galactic bulge , 
an ideal case for e xtinction studies (see also lAlves et aLl l2008l; 
lOnishi et al.lll999l) . We took this map as the true dust column 
density of a cloud complex, and simulated a set of background 
stars. In order to show the effects of substructures, we used a very 
small density of background stars (25 stars deg"^ instead of the 
original ~ 9.4 x 10^ stars deg"^ used to build the 2MASS/Nicer 
map). This, effectively, correspond to a linear downsize of the 
structures of the Pipe nebula by a factor ~ 60. 

The simulations were carried out using the following sim- 
ple technique. We generated a set of background stars uniformly 
distributed on the field of view. The stars were characterized by 
exponential number counts with exponent a = 0.34 in the K 
band and by intrinsic color and color scatters similar to the ones 
observed in the 2MASS catalog (see Table 2 of Lombardi 20051 
for a complete list of the parameters used). We added to each star 
magnitude the local value of the of the 2MASS/Nicer extinction 
and also some random measurement errors: 



/„ = J„ + kKAKiXn) + s\P , 

H„ = H„ + knAKiXn) + fif ^ 
K„ ^ K + A,(x„) + sf> , 



(36) 
(37) 
(38) 



where ej/'^''^^ denote random variables used to model the pho- 
tometric errors of the stars, and where as usual we used the 
hat to denote measured quantities. For simplicity, here, we took 
the errors as normal distributed random variables with constant 
variance 0.05 mag in all bands (the typical median variance of 
2MASS magnitudes); moreover, we assumed a hard complete- 
ness limit at 15 mag in all bands (i.e., if a magnitude exceeded 15 
we took the star as not detected in the coiTesponding band). We 
then applied the Nicer algorithm to these data, thus obtaining, 
for each star, its measured extinction and related error. Finally, 
we produced smooth extinction maps by interpolating the indi- 
vidual extinction measurements with three different techniques: 
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Fig. 5. Top left. The original 2MASS/Nicer map of the Pipe nebula used in the simulations. Top right. The mean reconstructed 
map that one would obtain if the density of background star were uniform. This is basically a smoothed version of the original 
map in the left. Levels are spaced at 0.2 mag. Bottom left. The difference between the average reconstructed maps, calculated from 
background stars following the density of Eq. and the average map at top right. Levels are spaced at 0.02 mag. Bottom right. 
Same map as bottom left, but using the improved method presented in this paper. Levels are spaced at 0.0025 mag. 
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Fig. 6. Left. Average reconstructed map using a 1-star (A^ = 1) Voronoi method, assuming a uniform distribution of stars. Levels are 
spaced at 0.2 mag. Right. Difference between the average reconstructed Voronoi map, obtained from background stars following 
the density of Eq. (|2]l, and the average true map to the left. Levels are spaced at 0.02 mag 




(1) simple moving average using Eq. ([T]i; (2) modified moving average maps A{x) obtained with the three techniques. We then 
average using Nicest [Eq. ([34]i1; (3) nearest neighbour interpo- compared these average maps with similar maps obtained from 
lation with a single star (A^ = 1). uniformly distributed stars (in order to produce such maps, 

we applied the completeness limit to un-extincted magnitudes). 
We repeated the whole procedure 1 000 times, using each j,^^^^^^ ^ ^ ^^^^ ^j^^ ^^^^^^^ obtained, which are in clear 
time a different set of random background stars, and took the 
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Fig. 7. Same as Fig.|5](bottom), but with a fraction of f - 0.05 of foreground stars. The increases of relative number of foreground 
stars in the denser regions of the cloud makes the bias more pronounced. Still, Nicest performs very well. Contours levels are spaced 
at 0.02 mag in the Nicer plot (left) and at 0.005 mag in the Nicest plot (right); contours at -0.1 mag or below are plotted in white. 



favour of Nicest. In particular, both the simple moving aver- 
age and the nearest neighbor suffer from significant biases, up to 
Ak ^0.1 mag or above, particularly in the most de nse regions 
of the Pipe nebula. As shown in va rious papers (e.g. lLada et alj 
119941; iLombardi et al.ll2006l l2008l) . the amount of small scale 
structure present in dark molecular clouds increases with the col- 
umn density, and thus it is not surprising that the larger biases 
are observed there. In contrast, the method presented in this pa- 
per reduces the bias by ~ 10, a factor comparable to the weight 
number of stars N (which in our simulations is N - 10 in the 
higher column density regions). 

The simulation performed here allows us to evaluate also 
the average quadratic difference between the extinction map ob- 
tained for the various methods A(x) and the true (smoothed) map 
A(x), i.e. the quantity 



{[A{x)-A(x)f) = Var[A(;«:)] + [{A{x) - Aix)x)f 



(39) 



As shown by the above equation, the quantity considered here 
can be written as the sum of the variance and the square of the 
bias of the extinction map. Our simulations show that, although 
Nicest, as expected, has a slightly larger variance then Nicer, 
there is still a gain by a factor at least two in the squared differ- 
ence of Eq. ( |39] |. In other words, the significant reduction in the 
bias performed by Nicest largely compensate the increase in the 
scatter introduced by this method. A naive comparison between 
the results obtained from Nicest and from the Voronoi method 
would provide even larger differences in favour of Nicest, but 
we stress that since the Voronoi method interpolates the extinc- 
tions over a fixed number of stars (typically smaller than the one 
employed by Nicest in our simulations) a direct comparison is 
not possible. 

Figure I?] shows the result of simulation carried out with the 
presence of foreground stars. In particular, we generated stars 
following the same prescriptions described in the previous para- 
graph, but allowing for a fraction /o = 0.05 of foreground ob- 
jects. Because of the effects of extinction, the effective fraction 
of foreground objects increases in the most dense regions of the 
cloud, which usually are also the ones with the larger substruc- 
tures: the two biases thus add up and can be particularly severe. 
As shown by Fig.|2l Nicest is effectively able to cope with rel- 
atively large fractions of foreground stars while still providing a 
virtually unbiased column density estimate. 
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Fig. 9. The difference between the modified extinction estimates 
[Eq. ( l28T ll and the standard ones [Eq. ([TJ] around Barnard 59. 
The contour levels are at A = [0.5, 1 .0, 1 .5, 2.0} mag of the map 
of Fig. E 



5. A sample application 

The method presented in this paper was finally applied to the 
whole 2MASS point source catalog for the Pipe nebula region. 
We used the 4.5 million stars of the 2MASS catalog located in 
the window 



-4° < Z < 4° 



+T <b< H-8° 



(40) 



The analysis was carrie d out following the prescriptions of 
'Lombar di & AlvesI (1200 1). but using the modified estimator ( l34l l 
to evaluate A. In particular, we generated the final extinction 
map, shown in Fig. [8] on a grid of about 1000 x 750 points, with 
scale 30 arcsec per pixel, and with Gaussian smoothing charac- 
terized by FWHM = 1 arcmin. The slope of the number counts 
was estimated to be or = 0.32 + 0.02 in the H band. 

The final, effective density of stars of about 8 stars per pixel 
guarantees that the approximation used to derive the unbiased es- 
timator ( [34l i is valid, and that a significant improvement over the 
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Fig. 8. The Nicest extinction map of the Pipe nebula, using the modified estimator (|28]) . 



0.30 




Fig. 10. The difference between the modified extinction esti- 
mates [Eq. (|28]|1 and the standard ones [Eq. ([T]i] around the 
peak ID 3 of iLombardi et alj (l2006l) . The contour levels are at 
A = {0.5, 1.0, 1-5, 2.0} mag of the map of Fig. [8l 



standard Nicer method can be expected. The largest extinction 
was measured close to Barnard 59, where A ^ 2. 68 mag (a value 
that is 0.41 mag larger than what obtained in (iLombardi et al.l 
2006)). 

Figures |9] and [To] show a comparison between dense regions 
mapped using Nicer and Nicest. We note that, as expected, the 
two methods are equivalent in low-density regions, while the 
new one consistently estimates larger column densities as A in- 
creases. The same effect can be appreciated more quantitatively 
from Fig.[TT] where we plot the relationship be tween the average 
estimates A obtained in ILombardi et al.l (l2006h and here. 

Final ly, we argue that the plo t of Fig.[TT]is strongly related to 
Fig. 9 of Lombardi et al.l (l2006h where we show the increase in 
the scatter of A„ for the different stars used in each point of the 
extinction map as a function of the average A in the point. This 
A-cr relation is most likely due to unresolved substructures in the 
cloud, that are expected to be more prominent in high-density 
regions. Recently, we re considered the A-cr rel ation and defined 
a quantity called A^(jc) (ILombardi et alj|2008h . This quantity is 
simply related to the local scatter of measured extinctions, but 
also properly takes into account the contribution of measurement 
errors to the observed scatter Interestingly, this quantity can be 
defined in terms of simple observables, and can be shown to be 
directly related to the local scatter of column densities: 



LnWn 



(41) 
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Fig. 11. The average difference AA between the estimators (l28T l 
and ([TJ as a function of the estimated column density A. The 
lower, solid line is relative to the Pipe nebula map at FWHM - 
1 arcmin resolution; the upper dashed line is at for a resolution 
FWHM = 1.5 arcmin. The latter plot is essentially what we 
would obtain at a FWHM = 1 arcmin if the distance of the Pipe 
nebula increases by a factor 1.5. 



where A„ = A{x,^ - A(x) is the difference between the column 
extinction at the position of the «-th star and the local average 
column extinction. A comparison of Eq. (1411 1 with Eq. dHJ clearly 
shows that, except for a numerical factor yS, the bias expected in 
the standard Nicer me thod is simply proportional to the l^{x) of 
jLombardi et al.ll2008l) . 

6. Discussion 

The numerical simulations and a first sample application of 
Nicest have shown that the method presented in this paper can 
significantly alleviate the bias introduced by small-scale struc- 
tures and by foreground stars in extinction studies. Of course. 
Nicest too has some limitations, which however are largely un- 
avoidable (and thus inherent to any extinction-based method). 

First, we note that in order for Nicest to work effectively, 
the weight number of background stars N must be significantly 
larger than unity: as shown in Sect. 13.21 N is directly related 
to the reduction in the bias provided by Nicest with respect to 
Nicer, and thus having N ~ \ does not give any benefit. This 
point is also important when correcting for foreground star con- 
tamination. For example, \i N ~ \ and the local fraction / 
of foreground stars is large (e.g., because we are in a particu- 
larly dense core), on average we do not expect any background 
star, and we will thus consistently measure a vanishing extinc- 
tion. Hence, Nicest, like any other extinction based method, only 
works if there are a sizable number of background stars that can 
be used for reliable extinction measurements. 

The above point might be interpreted as an exceedingly strin- 
gent requirement for the smoothing window used in Nicest. In 
fact, if a weight function 'w{x';x) - wix' - x") invariant upon 
translation is used, this function should be taken broad enough 
to guarantee that the weight number of background stars N{x) 
is significantly larger than unit everywhere. This, in turn, would 
imply that A/^ » 1 in the intermediate extinction regions, i.e. that 
the extinction map in these regions has a poor resolution, well 
below the limits imposed by the density of background stars. In 
reality, the whole derivation of the statistical properties of the 
Nicest technique is still valid for weight functions which are not 
spatially invariant. Hence, one does not need to use a fixed win- 



dow size for w(x'\ x), but rather this function could be taken to 
change shape for different locations x. For example, a simple 
scheme could be the use of a Gaussian shape for wix'^x) with 
the typical scale chosen according to a local estimate of the den- 
sity of background stars, in a way such that N(x) ~ const -10. 
This choice would guarantee an optimal resolution everywhere, 
and would still allow one to make use of the Nicer technique to 
reduce the inhomogeneity bias. 

Another point to keep in mind is that Nicest gives a large 
weight to red sources. This has two potential unwanted conse- 
quences. First, it introduces a small bias, that has been corrected 
to first order in Eq. ( |34l i. Second, as shown also by the numerical 
simulations described in Sect. |5] it slightly increases the noise in 
the final map. Again, since the noise in final map A(jc) is propor- 
tional to \ IN, the solution is to make sure that N is sufficiently 
large. We believe that a small increase in the noise is a fair price 
to pay for the significan t reduction in th e bias provided by Nicest 
over Nicer (see also .Eadie et al.|[T971h . 

A potentially more severe problem can be young stellar ob- 
jects (YSOs) with infrared excess. These sources can be present 
in the most dense star-forming cores of molecular clouds, where 
they can severely bias our method. We note, however, that since 
these sources have peculiar colors and are not usually present 
in the control field used for the calibration, they represent any- 
way a problem for all extinction based methods, including Nicer 
and the nearest neighbours ones. A median estimator is in prin- 
ciple safer to use in these cases, as long as the YSOs present in 
the core are a minority of the total number of background stars 
observed in the region; unfortunately, the tendency of YSOs to 
appear in clusters does not help here (note also that a median is 
usually more noisy than the simple average). In any case, the ob- 
vious solution is thus to exclude as much as possible YSOs from 
the list of sources used in the extinction maps. 

Finally, let us briefly mention alternative possibilities to re- 
duce the bias considered in this paper Since the substructure 
bias is due to a relationship between the local density p(jc) and 
the extinction map A{x') [Eq. (O], one could try to use a tracer of 
the local density to perform the correction. We discussed already 
a method based on this concept, the nearest neighbours interpo- 
lation. However, as shown in Sect. |4] this method fails (see also 
Sect. |2.2| i. In general, a problem with methods based on the local 
estimate of the density of background stars is that one needs to 
be able to obtain p{x) on scales smaller than the ones that char- 
acterize the weight w(x'\x), or otherwise it is not possible to 
give different weights to the different stars that contribute signif- 
icantly to the weighted average of Eq. ([T]l- However, if the den- 
sity of background stars is large enough to have a good handle on 
the local density p(x), one basically does not need to care about 
substructures, because it is always possible to make extinction 
maps at the same resolution as the density map. Hence, the only 
effective way to deal with substructures involves a use of the lo- 
cal extinction measurements, as done in Nicest. Note also that 
any density-based correction will be completely ineffective with 
foreground stars, in contrast to the method presented here. 

7. Conclusions 

The main results obtained in this paper can be summarized in the 
following points: 

- We discussed the effects of small scale inhomogeneities in 
the NIR extinction maps based on color excess methods, 
and showed that large inhomogeneities can significantly bias 
standard extinction maps toward low column densities. 
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- We proposed a new estimator for A, Nicest, and we showed 
that it is (i) equivalent to the usual estimator in low-density 
regions, (ii) unbiased and robust for any value of A, i.e. 
insensitive to the actual form and amount of substructure 
present in the cloud. 

- We showed that the new estimator is also suitable in presence 
contamination by foreground stars. 

- We tested Nicest against numerical simulations, and showed 
that it effectively reduces by a large factor the biases due 
to substructures and to foreground stars. We also tested an 
alternative approach, the nearest neighbor, and showed that 
the results obtained from this interpolation are still severely 
biased. 

- We applied Nicest to the Pipe nebula and showed a few pre- 
liminary properties of the resulting extinction map. We also 
noted a direct connecti on between the bias o f the Nicer and 
the map defined by jLombardi et alj|2008l) . 

Acknowledgements. We thank J. Alves, C.J. Lada, and G. Berlin for stimulating 
discussions and suggestions, and the referee, L. Cambresy, for helping us im- 
prove the paper This research has made use of the 2MASS archive, provided by 
NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion 
Laboratory, California Institute of Technology, under contract with the National 
Aeronautics and Space Administration. 



References 

Alves, J., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159 

Alves, J., Lombardi, M., & Lada, C. 2008, in Handbook of low mass star forma- 
tion regions, ed. B. Reipurt, ASP Conference Series 

Bok, B. J. 1937, The distribution of the stars in space (Chicago: University of 
Chicago Press, 1937) 

Cambresy, L., Beichman, C. A., Jan-ett, T. H., & Cutri, R. M. 2002, AJ, 123, 
2559 

Cambresy, L., Jarrett, T. H., & Beichman, C. A. 2005, A&A, 435, 131 
Cambresy, L., Petropoulou, V., Kontizas, M., & Kontizas, E. 2006, A&A, 445, 
999 

Dobashi, K., Bernard, J.-R, Hughes, A., et al. 2008, A&A, 484, 205 

Eadie, W., Drijard, D., James, E, Roos, M., & Sadoulet, B. 1971, Statistical 
Methods in Experimental Physics (Amsterdam New-York Oxford: North- 
Holland Publishing Company) 

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590 

Heyer, M. H. & Brunt, C. M. 2004, ApJ, 615, L45 

Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931 

Lada, C. J., Lada, E. A., Clemens, D. P, & Bally, J. 1994, ApJ, 429, 694 

Larson, R. B. 1981, MNRAS, 194, 809 

Lombardi, M. 2002, A&A, 395, 733 

Lombardi, M. 2005, A&A, 438, 169 

Lombardi, M. & Alves, J. 2001, A&A, 377, 1023 

Lombardi, M., Alves, J., & Lada, C. J. 2006, A&A, 454, 781 

Lombardi, M., Alves, J., & Lada, C. J. 2008, A&A in press 

Lombardi, M. & Schneider, R 2001, A&A, 373, 359 

Lombardi, M. & Schneider, R 2002, A&A, 392, 1 153 

Lombardi, M. & Schneider, R 2003, A&A, 407, 385 

Mathis, J. S. 1990, ARA&A, 28, 37 

Onishi, T., Kawamura, A., Abe, R., et al. 1999, PASJ, 51, 871 
Rieke, G. H. & Lebofsky, M. J. 1985, ApJ, 288, 618 
Vazquez-Semadeni, E. & Garcia, N. 2001, ApJ, 557, 727 
Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, ApJ, 161, L43 
Wolf M. 1923, Astronomische Nachrichten, 219, 109 



